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Abstract Theories of the gate-induced carrier density mod- 
ulation in bulk graphene are reviewed, and a numerical sim- 
ulation procedure based on electrostatic simulation using 
Matlab's pdetool is introduced. A brief introduction to the 
usage of the tool aiming at the present electrostatic prob- 
lem is introduced particularly for the readers who are as 
well Matlab users. The analytical capacitance models and 
the numerical simulation, namely the Poisson-Dirac itera- 
tion method, are compared, showing that the quantum ca- 
pacitance correction plays usually a minor role. Practical 
examples corresponding to realistic experimental conditions 
for graphene pnp junctions, superlattices, and linear back- 
ground potential are shown to illustrate the applicability of 
the introduced methods for calculating the carrier density 
modulation in bulk graphene. 

PACS 73.22.Pr • 85.30.De • 72.80. Vp • 41.20.Cv 



1 Introduction 

Electronic transport in graphene [1,2], a one-atom-thick 
honeycomb carbon lattice, is one of its main issues among 
the increasing number of fundamental studies ever since 
the first successful isolation of stable monolayer graphene 
flakes in 2004 [3]. What led to the explosive growth of 
the graphene literature, however, was not only the discov- 
ery of the mechanical exfoliation (Scotch-tape method) for 
graphene flake preparation, which made graphene easily 
accessible to laboratories all over the world, but also the 
characterization of the electronic properties of graphene by 
electrical gating, which provides a direct way to modulate 

Ming-Hao Liu (IPJJE) 

Institut fur Theoretische Physik, Universitat Regensburg, 
D-93040 Regensburg, Germany 
E-mail: minghao.liu.taiwan@gmail.com 



the carrier density, and hence the Fermi level, of graphene 
[3]. Conductance (resistance) sweep using a single backgate 
is therefore a standard electronic characterization tool for 
graphene. Double-gated graphene opens even more possibil- 
ities of graphene electronics and allows experimental stud- 
ies of graphene pn and pnp junctions [4,5,6,7,8], as well as 
the interesting physics of Klein tunneling [9, 10, 1 1, 12, 13]. 
Gate-induced carrier density modulation therefore plays an 
essential role for fundamental as well as advanced studies of 
graphene electronics. 

Theory of the gate-induced carrier density modulation is 
mainly an electrostatic problem. How one should obtain the 
gate-voltage dependence of the carrier density in graphene 
depends actually on how precise one wishes. For cheapest 
computation, the graphene sheet carrier density can be di- 
rectly regarded as the induced surface charge density adja- 
cent to graphene [3], which is treated as a conductor fixed at 
zero potential. This corresponds to the classical capacitance 
model (CCM) that is widely adopted in most experimental 
works on graphene transport [2] and can be solved exactly. A 
more precise computation takes into account the relation be- 
tween the induced charge density on graphene and the elec- 
tric potential energy that those charge carriers gain, through 
the graphene density of states [14,15,16]. This makes the 
computation iterative [8, 17, 18] and hence a bit more expen- 
sive, but actually corresponds to the quantum capacitance 
model (QCM) [19], where an exact solution for single-gated 
graphene is possible [16]. Further consideration such as the 
Coulomb interaction of the induced charges on the graphene 
sheet is possible [15, 17], but this would be out of the scope 
of the present discussion. 

Whereas a thorough and comprehensive review on the 
theory of gate-induced carrier density modulation of bulk 
graphene so far does not exist in the literature, the major 
part of this paper aims at providing this missing piece. The 
review includes both the analytical and numerical aspects, as 
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well as a brief introduction to the usage of Matlab's pdetool, 
in order for a self-contained context. Readers who happen 
to be Matlab users would find this brief usage helpful, but 
non-Matlab users may as well neglect it without encounter- 
ing further gaps. The analytics based on the CCM and QCM 
and the numerics based on the iterative method, namely, the 
Poisson-Dirac method (PDM), using Matlab's pdetool will 
be compared, showing that the quantum correction plays 
usually a minor role, unless the oxide between graphene and 
the metallic gate is exceedingly thin. In the case of single- 
gated graphene, consistency between the QCM and PDM 
is satisfactory even for nonuniform capacitors and becomes 
exact for uniform capacitors. 

With a full understanding of the gate-voltage modula- 
tion of the graphene carrier density, examples of its appli- 
cations aiming at providing realsitic local energy band off- 
sets due to electric gating will be discussed. This is partic- 
ularly important for an accurate computational modeling of 
transport in graphene. Examples include the cases of (i) a 
single graphene pnp junction, (ii) a graphene superlattice, 
and (iii) a background potential linear in position. Practi- 
cally, the case of (i) provides the study of the physics of 
Klein backscattering [20,12,21], while the combination of 
the cases (ii) and (iii) is the underlying prerequisite of the 
Bloch-Zener oscillation in graphene [22]. 

This paper is organized as follows. In Sec. 2, we first 
provide a brief introduction to the usage of Matlab's pdetool, 
pointing out the least knowledge required to apply the tool 
on the present specific electrostatic problem. Theories of the 
gate-induced carrier density modulation in graphene are re- 
viewed in Sec. 3, where the analytics based on the capac- 
itance models and the numerics based on the PDM is also 
compared. Practical applications based on the theories re- 
viewed in Sec. 3 are given in Sec. 4, and a summary of the 
present work is concluded in Sec. 5. 

2 Usage of Matlab's pdetool for electrostatics 

The pdetool is a useful numerical tool built in Matlab and 
provides a convenient way to solve several classic partial 
differential equation (PDE) problems in two-dimension. For 
the electrostatics at our present interest, the Poisson equa- 
tion, 

-V.(e r V M ) = -, (1) 
£o 

obtained from two of the Maxwell's equations, V x E = 
and V - D = p, which respectively lead to E = — Vu and 
V • £, E = p/Eo, is the central equation that the pdetool solves 
for the electric potential u} In Eq. (1), the product of the di- 

To be consistent with the pdetool, we name the electric potential as 
u, while reserve the variable V for the energy band offset (the "on-site 
energy" in the language of tight-binding formulation). 
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electric constant (relative permittivity) e r and the free space 
permittivity £o gives the absolute permittivity e = £,£o. 

A full introduction to the usage of the pdetool can be 
found in the Matlab documentation [23] and need not be re- 
peated here. To digest the full user's guide of the tool, how- 
ever, is not necessary for our focussed problem, which is es- 
sentially an electrostatic problem. This section is basically 
to elaborate those that are less clear in [23] but nevertheless 
important for our purpose of obtaining the gate-voltage de- 
pendence of the graphene carrier density, and to point out 
the least required knowledge for this purpose. 

2. 1 Overview of pdetool 

To solve a PDE problem using the pdetool, required neces- 
sary inputs can be exported from the graphical user interface 
(GUI) of the pdetool (initiated by executing "pdetool" from 
the command window) and are briefly described in the fol- 
lowing. 

(i) System geoemtry. The geometrical shapes of the build- 
ing blocks, such as the oxide layers, metallic gates, etc., 
which constitute the system where the PDE problem is 
defined, can be drawn in the "Draw Mode" of the GUI. 
The resulted "decomposed geometry" allows us to pro- 
ceed to the rest of the inputs, but there is no need to 
"Export Decomposed Geometry, Set Formula, Labels..." 
from the "Draw menu" since they will not be needed by 
the PDE solvers. 

(ii) PDE coefficients. In the "PDE Mode" of the GUI, one 
can designate different regions of materials by filling in 
the respective dielectric constants and space charge den- 
sities. These are stored in certain PDE coefficients ma- 
trices, which can be output from the GUI and will be 
required by the PDE solvers in programming. 

(iii) Boundary conditions. In the "Boundary Mode" of the 
GUI, boundary conditions for each bounding edge can 
be assigened. The resulting boundary matrix b, which 
will be required by the PDE solvers in programming, 
and the Decomposed Geometry, which will be required 
when visualizing the PDE geometry, can be exported 
from the "Boundary menu". An elaborated instruction 
about b will be given later. 

(iv) Mesh points. The mesh points are those spatial points 
at which the numerical solutions are desired. They can 
be created, refined, or jiggled in the "Mesh Mode" of 
the GUI. The resulting triangular mesh data, stored by 
point, edge, and triangle matrices, can be exported by 
the GUI and will be used not only when calling for the 
PDE solvers but also when visualizing the solution. 

With all these requirements completed, the PDE prob- 
lem is then defined, and the solution can as well be ob- 
tained by clicking "Solve PDE" within the GUI, which is 
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user-friendly but cannot be "programmed". When perform- 
ing certain real calculations, however, especially when a sys- 
tematic change of variables is required, programming with, 
e.g., looping, is inevitable and the requirements of (ii)-(iv) 
will be the necessary inputs of the PDE solvers. For our pur- 
pose of simulating the carrier density modulation due to gat- 
ing, we would often need to change the gate voltages, which 
are described by the boundary conditions. Thus although 
each of (ii)-(iv) can be programmed by using relevant com- 
mands, in the following only the implementation of (iii) by 
commands will be described in detail. 

2.2 Boundary conditions 

2.2.7 The boundary condition matrix 

By searching "assemb" from the Matlab help, we see that 
the boundary conditions are saved in a matrix called b, with 
the following data format: 

• Row 1 contains the dimension N of the system. (Note: 
normally N = 1. If one solves two coupled variables, 
then N = 2, etc.; by examining the exported boundary 
condition matrix b, one woule find that N = for inner 
boundaries.) 

• Row 2 contains the number M of Dirichlet boundary 
conditions. 

• Rows 3 to 3 + N 2 — 1 contain the lengths for the strings 
representing q. The lengths are stored in columnwise or- 
der with respect to q. 

• Rows 3 + N 2 to 3 + N 2 + TV - 1 contain the lengths for 
the strings representing g. 

• Rows 3 + N 2 +N to 3 + N 2 +N + MN - 1 contain the 
lengths for the strings representing h. The lengths are 
stored in columnwise order with respect to h. 

• Rows 3+N 2 +N + NM to 3 + N 2 +N + MN + M - 1 
contain the lengths for the strings representing r. 

• The following rows contain text expressions represent- 
ing the actual boundary condition functions. 

Here, two types of boundary conditions 2 are included, 
namely, the Neumann boundary 

n- (cVm) +qu = g, (2) 
and the Dirichlet boundary 

hu = r. (3) 

Here, c contains the PDE coefficients (such as the dielec- 
tric constant in different regions), and n is the normal of the 
boundary. So the boundary conditions for given gate volt- 
ages would be the Dirichlet type, with h = 1 and r being 

2 The mixed type boundary conditions will not be encountered in 
the present discussion. 



the corresponding voltage. For the Neumann type boundary 
condition, we normally consider q = 0, and g represents the 
surface charge. 

In the following, let us be more specific about the format 
of the boundary conditions matrix, considering the two types 
of boundaries with N = 1 . 

2.2.2 Dirichlet boundary 

The format of the boundary matrix b is described as follows. 

• Row 1 contains the dimension TV of the system: 1 . 

• Row 2 contains the number M of Dirichlet boundary 
conditions: 1. 

• Row 3 contains the length for the strings representing q, 
which is 1 since q = 0, though not used. 

• Row 4 contains the length for the strings representing g, 
which is 1 since g = 0, though not used. 

• Row 5 contains the length for the strings representing h, 
which is 1 since h=l. 

• Row 6 contains the length for the strings representing r. 

• Then comes the text expressions of q,g,h,r. 

An example of a Dirichlet boundary with r — 3.5 would 

be: 

b = [1 1 1 1 1 3 '0' '0' '1> '3.5'] ' ; 

The boundary condition may include the x and y position 
coordinates and their functions, and can be written even in 
terms of the solution u (nonlinear solver required). For ex- 
ample, 

b = [1 1 1 1 1 4 '0' '0' '1' 'x.~2'] ' ; 
For another example, 

b = [1 1 1 1 1 9 '0' '0' '1> 'sin(x) .~u>] ; 

2.2.3 Neumann boundary 

The format of the boundary matrix b is described as follows. 

• Row 1 contains the dimension TV of the system: 1 . 

• Row 2 contains the number M of Dirichlet boundary 
conditions: 0. 

• Row 3 contains the length for the strings representing q, 
which is 1 since q — 0. 

• Row 4 contains the length for the strings representing g. 

• Then comes the text expressions of q,g. 

An example of a Neumann boundary with surface charge 
g = 1 .6 would be 

b = [1 1 3 '0' '1.6'] ' ; 

Another example 

b = [1 1 21 '0' '-13.295*sign(u) .*u. "2'] ' ; 

will actually be used when we apply the Poisson-Dirac iter- 
ation method. 
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2.2.4 Text expression of the boundary conditions 

The boundary condition matrix b exported from the GUI 
of the pdetool looks filled with purely numbers. This is the 
"number representation" of the text strings. For example, a 
number 48 within the b matrix actually means '0' : 

» char (48) 

ans = 



Conversely, if we want to transform the strings into num- 
bers, we can simply use the 'double' command: 

>> double ( 'x. "2 ') 
ans = 

120 46 94 50 

Thus to enter a boundary condition of, e.g., — (x 2 + y 2 ), we 
may fill in: 

double ('-(x.~2+y.~2) ') ' 

To enter a Dirichlet boundary condition of a given number 
assigned by a variable named Vtg, we may fill in: 

double (num2str (Vtg) ) ' 

Note that in the above two examples, the operator ' at the 
end is to take transpose of the converted text strings since 
the boundary conditions are saved columnwise in b. 

For a real PDE problem, the number of columns of the b 
matrix depends on the total number of edges, including inner 
and outer boundaries. The «th column records the boundary 
condition for the «th edge. Thus before exporting the initial 
b matrix from the GUI of pdetool, one has to check the cor- 
responding boundary label (by showing the edge labels in 
the "Boundary Mode"). 

2.3 Some important commands 
2.3.1 Solving the PDE 

A standard PDE solver is called assempde. An example for 
its usage is as follows. 

u = assempde(b,p,e,t,c,a,f ) ; 

°/ b: matrix of boundary conditions 

°/ p: points of the mesh grid 

7. e: edges 

7« t : triangles 

% c,a,f: coefficients of the pde problem 



When the solution is involved in the boundary condi- 
tions, the solution mode has to be switched to nonlinear. An 
example for its usage is as follows. 

[u.res] = pdenonlin(b,p, e , t , c ,a,f , ... 

'report' , 'on' , 'Maxlter' ,le5, 'u0' ,u0) ; 
7, u: the solution, res: not important here 
7o p,e,t,c,a,f : same as above 
7o 'Maxlter': maximal number of iter, rounds 
7o 'u0': initial guess of the solution 

2.3.2 Interpolation 

To find the values at those points one desires, an important 
command called tri2grid should be used, which interpo- 
lates from the PDE triangular mesh to a given rectangular 
grid. An example for usage is as follows. 

uxy = tri2grid(p,t,u,x,y) ; 
7o u: the obtained solution 
7. p,t: same as above 

7» x,y: rectangular grid points for interpol. 
2.4 Remarks on units 

In the pdetool, everything is displayed with dimensionless 
numbers. The actual units can be chosen as what we would 
like. Deducing relevant coefficients for a specific set of cho- 
sen units is therefore important before we use the pdetool to 
solve any actual problems. 

In addition to the physical units, the free space per- 
mittivity £o is suppressed throughout the program. Recall 
the boundary conditions for the displacement field D at a 
conductor-dielectric boundary, which can be derived by ap- 
plying the Gauss's law: D t = and D„ = p s , where D t and 
D n represent the components tangential and normal to the 
interface, respectively. The normal component of the dis- 
placement field D n therefore means the surface charge den- 
sity: 

-e r V M -n=— . (4) 



Comparing Eq. (4) with Eq. (2) with q = and c represent- 
ing e, (instead of e), one can see that the "Surface charge" 
actually means 




when filling the Neumann boundary condition in the 
"Boundary Mode" of the GUI. Similarly, when filling the 
PDE coefficients in the "PDE Mode", the "Space charge 
density" actually means rho = p/eo, i.e., the right-hand 
side of Eq. (1). 
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3 Theories of gate-induced carrier density modulation 

In this section, analytical theories of the gate-induced car- 
rier density modulation, including the classical and quan- 
tum capacitance models, are briefly reviewed, a numerical 
scheme of the Poisson-Dirac iteration method is introduced, 
and a numerical comparison between analytics and numer- 
ics is provided. 



3.1 Classical capacitance model 
3.1.1 The model 

We begin with the classical capacitance model, which con- 
siders a parallel-plate capacitor composed of an oxide di- 
electric with permittivity e sandwiched by a metallic gate 
(at z — d) and a conducting graphene sheet (at z — 0) as 
sketched in Fig. 1(a). Let the electric potential at the gate 
be u(x,z = d) — V g and the graphene layer be grounded: 
u(x,z = 0) = Vq = 0. The surface charge density at z — 
+ (the surface of the oxide dielectric in contact with the 
graphene layer) from Eq. (4) is given by 



du 

P.s = -«J 
OZ 



e d -Q- LoxV >' 



(6) 



where C ox = B/d is the classical capacitance (per unit area) 
of a uniform parallel-plate capacitor. Regarding this surface 
charge (6) directly as those in the graphene layer, we have 
the carrier density 



P.S Cos 

nc = — = — Vg 



(7) 



which is a widely used formula for estimating the graphene 
carrier density [2]. For uniform capacitors with C ox = B/d, 
Eq. (7) numerically reads 



BV g B r V 



n c = —r = x 5.5263 x 10 12 ctrr z , 
ed d 

where V g and d are in units of V and nm, respectively. 

3.1.2 Using pdetool 

The Dirichlet boundary conditions 



(8) 



u(x,z) 



0, at graphene boundary 
V g , at gate boundary 



(9) 



can be straightforwardly implemented. 3 The standard PDE 
solver assempde introduced in Sec. 2.3 should be chosen. 

3 If a uniform capacitor (without x dependence) is desired, one needs 
to assign Neumann boundary conditions at the left and right sides of the 
oxide boundaries with vanishing surface charge density g = 0, which 
forces the displacement field to be normal (tangential) to the side (top 
and bottom) boundaries. 



gate (g) 



\V = V a 



+ + + + + + + + + + + 
p = +en 



oxide (ox) 



z = 



p = —en 



graphene (G) 

(a) 

V 



V = V G 



Co 



V G 




(b) 



(O 



Fig. 1 (a) Schematic of a single-gated graphene. (b) Equivalent circuit 
plot of the classical capacitance model, (c) Equivalent circuit plot of 
the quantum capacitance model. 



Working with units V and nm, the carrier density (7) is nu- 
merically given by 



n c (x) = £, 



du(x,z) 



dz 



x 5.5263 x 10 12 cm -2 . 



(10) 



Note that the interpolation command tri2grid introduced 
in Sec. 2.3 may be useful in performing the numerical 
derivative du/dz at z = 0. 

3.1.3 Remark on the gate-induced Rashba spin splitting 

At this stage we may also estimate for graphene the gate- 
induced Rashba spin splitting, an intrinsic coupling between 
the spin and orbital degrees of freedom of charge carriers 
in a two-dimensional conducting plane subject to a per- 
pendicular electric field [24,25]. In graphene, the Rashba 
spin splitting has been shown by first principles to exhibit 
a linear dependence on the electric field [26,27]: Ar 
0.01 |J5|meV, where E is the electric field strength perpen- 
dicular to graphene given in units of V / nm. If the graphene 
carrier density n stems from gating, the corresponding sur- 
face charge density \p s \ = e\n\ = E \E\, in fact, has already 
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revealed the displacement field on itself, allowing us to ex- 
press the Rashba spin splitting in terms of the carrier density, 



A R = — x 1.8095 x l(T eV, 



(11) 



where n is in units of 10 12 cm~ 2 , a typical order of the 
graphene carrier density. 

This estimation indicates that the Rashba spin splitting 
induced solely by electric gating typically lies in the order 
of /leV, which may hinder the observation of those inter- 
esting physics based on the Rashba spin-orbit coupling in 
graphene, such as the interfacial spin and charge currents 
[28,29], or the spin-dependent Klein tunneling [30,31]. The 
position dependence of the Rashba coupling across a pn 
junction interface [32], on the other hand, can be accurately 
taken into account by putting the jc-dependence of n (or even 
£ r ) in Eq. (1 1). A stronger Rashba spin splitting in graphene 
is therefore less possible by gating, but may be achieved by, 
for example, using a ferromagnetic substrate with an inter- 
calated gold monolayer [33]. 



linear response regime, gains an energy — eVc from the elec- 
tric field, where — e is the electron charge 4 and Vq is the elec- 
tric potential at the graphene sheet obtained by solving the 
Poisson Eq. (1). This potential energy — eVc, which is equiv- 
alent to the "on-site energy" in the tight-binding transport 
formulation (see, for example, [21]), will raise the whole 
band structure, and thus lower the Fermi level by the same 
amount. We can therefore legitimately put 



E = -(-eV G ) = +eV G 
into Eq. (12), leading to 

P.v —en 



= ~ — s g n ( V g)z7Z — XT- 
e K(hv F ) 1 



(13) 



(14) 



The surface charge density at the graphene layer is now ex- 
pressed in terms of the solution u(x,z = 0), but is at the same 
time the Neumann boundary condition that influences the 
numerical solution to the Poisson equation. This formally 
makes the solution process iterative. 

3.2.2 Using pdetool 



3.2 Poisson-Dirac iteration method 

From Eq. (6) to Eq. (7), the assumption that "the in- 
duced surface charge density at the dielectric surface is the 
graphene carrier density" obviously have neglected a few 
physical details, such as the graphene density of states that 
govern the statistics of how the states in graphene should be 
filled by the carriers accordingly. In addition, filling the car- 
riers into graphene causes the change of its Fermi level, im- 
plying a potential energy shift that should further correspond 
to the electric potential times the electron charge. These are 
what the classical capacitance model have neglected and 
what the following Poisson-Dirac iteration method is going 
to compensate. 

3.2.1 Basic idea 

Consider a pristine graphene with Fermi level lying exactly 
at the charge neutrality point, i.e., the Dirac point Ep = 0. 
Application of the gate voltage V g induces a certain amount 
of additional charges on graphene, p v = —en, which oc- 
cupy the states in graphene according to its density of states 
D{E) = 2\E\/n (hv F f (within the Dirac model): 



n(E)= / D(E')dE' = sgn(£) 



% \hv F ) 



(12) 



where vp « 10 8 cm/ s is the Fermi velocity in graphene. A 
positive (negative) electron number density n raises (lowers) 
the Fermi level from to E. On the other hand, the electron 
at the Fermi level, which is responsible for transport in the 



The Dirichlet boundary condition (9) for the gate boundary 
remains valid, while that for the graphene boundary has to 
be modified to the Neumann type: 



u(x,z) 



at graphene boundary 
at gate boundary 



(15) 



where g = p s / £o is given by Eq. (14). Working with units V 
and nm together with vp = 10 8 cm / s, Eq. (14) becomes 



£ = -13.295sgn(V G )(^) — , 
£o V V J nm 



(16) 



which should be keyed as "-13 . 295*sign(u) . *u.~2"in 
the boundary condition matrix, noting that the solution in 
the pdetool is by default named u. The nonlinear solver 
pdenonlin introduced in Sec. 2.3 has to be chosen in this 
case, where the solution is involved in the boundary condi- 
tions, and the iteration will be automatically processed by 
the pdetool. 

Once the solution u(x,z), and hence the electrostatic po- 
tential at the graphene layer Vq(x) = u(x,z = 0), is itera- 
tively obtained, the desired carrier density profile n(x) can 
then be expressed in terms of Vg(x): 



n PD {x) =7.3471 x 10 13 x sgn[V G (Y)] 



Vg(x) 
V 



(17) 



which follows from Eqs. (12) and (13). Note that we have 
added explicitly a subscript PD in Eq. (17) to distinguish 
with the classical contribution, nc. 



Throughout this paper, e - 
tary positive charge. 



1.60217733 x 1(T 19 C is the elemen 
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3.3 Quantum capacitance model 



3.3.2 Quantum correction for uniform capacitors 



The relation between the induced charge density on 
graphene and the electric potential energy that those charge 
carriers gain through the graphene density of states is taken 
into account in the PDM, with a price of iteration process 
paid. For single-gated graphene, there is an alternative that 
can take this into account analytically: the quantum capac- 
itance model [19], which we briefly review here for bulk 
graphene following the work of [16]. 

3.3.1 The model 

The single-gated graphene shown in Fig. 1(a) is treated by 
the equivalent circuit plot as shown in Fig. 1(c), where an 
additonal capacitor Cq is inserted between the voltage point 
Vg and the ground, as compared with the CCM, Fig. 1(b). 
As Fig. 1(c) suggests V g = V G + V ox , using C ox = \p\ /V ox = 
en /Vox we have 



v K = v G - 



Co 



Co 



classical 



-—V G 



quantum 



(18) 



Following the same physics stated in Sec. 3.2.1, the carrier 
density at the graphene layer, i.e, Eq. (14) divided by — e, is 
expressed in terms of the electric potential thereof as 



n=sgn(eyc) K£) 



(19) 



Equating (18) and (19), one obtains a quadratic equation for 

Vg, 



sgn(eVc) 



■Vg\ 2 = Co 



K \h\>f I 



■v g -- 



Co 



-V G . 



Solving Eq. (20) for Vg and putting back to Eq. (18), the 
graphene carrier density can be written as 

n = tic + An, (21) 

where nc given by Eq. (7) is the classical contribution, and 



An 



n Q ( 1 -sgn(nc) J 1+2-^- 



with 



n ( C 0X hv F 



(22) 



(23) 



corresponds to the quantum correction. 5 

5 Note that Eqs. (21)— (23) were first derived in [16] and reviewed in 
[2], but a factor of 2 in the square root of Eq. (22) corresponding to the 
formula (1.15) in [2] was missing. 



For uniform capacitors, the well-known oxide capacitance is 
given by C ox = £ /d, or, by definition written as 



C 



e-n c 



(24) 



such that Eq. (23) and hence the quantum correction (22) 
is solely determined by the classical contribution nc- In this 
case we can further express Eq. (23) as ng = (e r /d) 2 x 2. 
0784 x 10 n cm~ 2 , where d is in units of nm and \>f = 
10 8 cm / s is again adopted. Together with Eq. (8), the quan- 
tum correction given by Eq. (22) can be written as 



An = sgn(V g )(^Y 



V„ d 



(25) 



x 2.0784 x 10 n cirT 2 



where V g is in units of V. We will soon see that this analyti- 
cal QPM for a uniform capacitor will exactly correspond to 
the numerical PDM. With nonzero V g and large d, one can 
further approximate An as 



An k, -sgn(y^) 



V g \ x 1.515< x 10 12 ctrT 2 



(26) 



which implies a rapid decay of An with d. 



(20) 3.3.3 Remark on quantum capacitance 



Note that the appearance of Cq stems from the finite density 
of states provided by the conducting layer for the electrons 
to occupy following the quantum nature of the Pauli ex- 
clusion principle, and hence the name quantum capacitance 
[19], which is not restricted to the material graphene. The 
expression of Cq for graphene [16], however, is not impor- 
tant for the present discussion. Instead, Cq leads to a quan- 
tum correction to the gate-induced carrier density An, which 
is the main focus here. 

In addition, recent experimental progress on the mea- 
surement of graphene quantum capacitance [34, 35, 36] sug- 
gests that the electron-hole puddles [37] induced by charged 
impurities may influence Cq at energies close to the charge 
neutrality point. The corresponding carrier density fluctua- 
tion Sn, which can be considered to develop a microscopic 
model to account for the smoothing of the graphene quan- 
tum capacitance at the charge neutrality point [38], is be- 
yond the scope of the present discussion. 
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3.4 Analytics vs numerics 

The two analytical capacitance models and the numerical 
scheme of the Poisson-Dirac iteration method are compared 
in the following, considering a single-gated graphene with 
individually infinite and finite sizes of the gate. 

We start with a uniform capacitor with infinitely ex- 
tending parallel conductors with different spacings d = 
100,200,300 nm. The oxide dielectric is regarded as S1O2 
with £,. = 3.9. The classical contribution nc is first com- 
puted following Sec. 3. 1 .2 [which gives the same result with 
Eq. (8)], and the quantum correction An is computed in two 
ways. For the analytical QCM, Eq. (25) is used to compute 
An. For the numerical PDM, the full carrier density «pd is 
computed following Sec. 3.2.2, and the correction is given 
by the difference «pd — «c- As shown in Fig. 2, the corre- 
spondence between the two approaches is exact. 

Next we consider a finite-size suspended topgate (such 
as those fabricated in [7]) with various voltages V tg = 
5,10,15 V and an extremely narrow spacing d = 20 nm. 
The calculations are similar to those for the uniform case 
described above. The only difference is the approximating 
form of the dielectric capacitance, 

C 0X M = ^, (27) 

from which hq given by Eq. (23) and hence the quantum 
correction Eq. (22) from the analytical QCM is obtained. 
The total carrier densities calculated by the iterative PDM, 
rcpo, and that by the QCM, hqc, are compared in the left 
panel of Fig. 3; the difference can hardly be told. In the right 
panel of Fig. 3, the total carrier density is compared with the 
classical contribution ric\ the difference is merely moderate, 
implying that «c lS still the main contribution to the gate- 
induced carrier density. The quantum corrections from the 




-200 -100 100 200 -200 -100 100 200 

x (nm) x (nm) 



Fig. 3 Left: Comparison between the position-dependent total carrier 
density obtained from the Poisson-Dirac method rcprj(jr) and that from 
the quantum capacitance model ;!qc(x), at various topgate voltages. 
Right: Total carrier density /ipoM vs the classical contribution nc(x). 
Inset: Spatial distribution of the electric potential V(x,z) due to the 
presence of the topgate with cross section 50 nm x50 nm suspended 
20 nm above the graphene sheet. The length unit is nm. 

two approaches are compared in Fig. 4, where a clear but 
still moderate discrepancy can be seen. 

From the above testing calculations, we may conclude 
that the classical contribution uq always plays the dominant 
role, and the QCM is equivalent to the PDM in the case of 
uniform capacitors. For nonuniform capacitors, the approx- 
imation of Eq. (27) that is generalized from Eq. (24) allows 
the non-iterative QCM to work equally well as compared to 
the iterative PDM, for the case of single- and finite-gated 
graphene. 

3.5 Beyond single gating 

The above discussion considers only a single-gated 
graphene. For double-gated graphene with topgate and back- 
gate at two sides of the graphene layer, the two gates can 
be regarded as independent, and their contributions to the 
carrier density modulation can be treated separately. When 
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Fig. 2 The quantum correction to the gate-induced earner density on 
graphene, calculated by the analytical quantum capacitance model and 
the numerical Poisson-Dirac iteration method. Inset: Two-dimensional 
plot for the iterative solution u(x,z) with length unit nm. 
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Fig. 4 Comparison of the correction to the carrier density from the 
Poisson-Dirac method, npo(x) — nc(x), and that from the quantum ca- 
pacitance model, An(x) from Eq. (22), at various topgate voltages. 
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multiple gates are fabricated at the same side, however, such 
as using an embedded local gate in addition to a global back- 
gate to creat ballistic pnp junction [13] (see Sec. 4.2), or 
patterned topgates that may generate a graphene superlattice 
(see Sec. 4.3), they should be simultaneously treated. 

In fact, the CCM (Sec. 3.1.2) as well as the PDM (Sec. 
3.2.2) are not restricted to the case of single-gated graphene. 
These two approaches work for any kind of gating geome- 
try, provided that the boundary conditions [Eq. (9) for CCM 
and Eq. (15) for PDM] at the graphene layer are properly 
assigned. The applicability of the QCM (at the level that we 
introduced above), however, depends then on the gating ge- 
ometry. In certain cases where Eq. (27) becomes ill-defined, 
the model cannot be simply applied. Generalization of the 
model to take into account composite gating geometry is 
possible, but is beyond the scope of the present discussion. 



4 Applications 

A successful simulation for electronic transport in bulk 
graphene relies on not only sophisticated computation tech- 
niques but also a realistic "potential profile" V that is experi- 
mentally relevant [21]. The term "potential" refers to the lo- 
cal potential energy added to the system Hamiltonian when 
modeling for graphene electronic transport. Thus the poten- 
tial profile simply means the local energy band offset of the 
graphene sheet subject to a spatially varying carrier density 
due to electrical gating. This section is devoted to the appli- 
cation of the carrier density calculation: the corresponding 
potential profile, or the local energy band offset, which is 
a simple computational task but nevertheless important for 
graphene electronic transport calculations. A few concrete 
examples will be illustrated, after a short review of the po- 
tential profile is given. 



4. 1 Potential profile (local energy band offset) 

In Sec. 3 we have introduced how to compute with a sat- 
isfactory accuracy the graphene carrier density n, which 
is related to the quasi-Fermi level through Eq. (12) as 
Ef = sga(ii)hvf y/n\n\. The energy band offset then reads 
V = E F — Ef, where E F is the global Fermi level. Choos- 
ing E F = (as is usually the case and will be adopted in 
the rest of the calculations), the space-resolved band offset 
reads [21] 



graphene 



V(x) = —sgn[n(x)]hvfy/it\n(x)\ 



= - U . 667 x sgn[«(*)] y meV 



(28) 
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which is termed on-site energy in the tight-binding formula- 
tion for transport calculations. 



Fig. 5 A graphene pnp junction using a global backgate with various 
voltages Vb e and an embedded local gate fixed at V\ g = 4 V. (a) The it- 
erative solution u(x,z) in units of V to the electrostatic potential within 
the oxide, subject to Vj,g = 30 V. The induced graphene carrier densities 
(b) and the corresponding potential profiles (c) based on the classical 
capacitance model and the Poisson-Dirac method. 



It should be remarked that the Poisson-Dirac iterative so- 
lution to the electric potential at the graphene layer Vq times 
— e readily gives the desired energy band offset, and Eq. (28) 
is not needed within this approach. Within the capacitance 
models, however, the electric potential on graphene is fixed 
at zero as a Dirichlet boundary condition, and Eq. (28) is 
needed. In other cases where PDM is partly used but the 
total carrier density is separately computed, one needs Eq. 
(28) as well. 



4.2 Graphene pnp junctions 

We begin the illustrative examples with a graphene pnp junc- 
tion, using a global backgate and an embedded local gate. 
The gating geometry is sketched in Fig. 5(a), simular to 
those experimentally fabricated in [13]. In this case both 
of the global and local gates influence the graphene carrier 
density from the same side, and therefore have to be treated 
at the same time. The QCM does not apply here since Eq. 
(27) is ill-defined, but nevertheless can be used to estimate 
the quantum correction due to the embedded local gate at 
the region close to it. We will mainly compare the results 
from the CCM and those from the PDM, fixing the local 
gate voltage at Vj g = 4 V while varying the backgate voltage 
with V hg = -60, -30, 0,30 V. 

The computed carrier densities are shown in Fig. 5(b), 
where the Poisson-Dirac solution agrees well with the sim- 
ulation presented in [13]. Since the local gate is embed- 
ded only lOnm below the graphene sheet, the quantum cor- 
rection excluded in the CCM becomes pronounced within 
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the locally gated region, and can be estimated by Eq. (25) 
of the QCM. At the center of the locally gated region, 
the electric field generated by the backgate is almost com- 
pletely screened, and the classical contribution to the car- 
rier density can be estimated by nc{x = 0) = eV\g/ed = 8. 
6211 x 10 12 cm~ 2 , leading to n e = 3. 161 2 x 10 I0 cm- 2 and 
hence An = —7.0735 x 10 n cm~ 2 , which is prettey close 
to « PD (x = 0) - n c (x = 0) = -7.0746 x 10 1 1 cm~ 2 from the 
data of Fig. 5(b) for all V hg . 

The carrier density profiles n(x) of Fig. 5(b) are trans- 
lated into V(x) via Eq. (28), as shown in Fig. 5(c). The pos- 
itive Vi g charges the locally gated graphene with a positive 
number of electrons, forming an n-type region with positive 
quasi-Fermi level Ef(x) > that is equivalent to applying 
a negative energy band offset V(x) < 0. Outside the locally 
gated region, the carrier type of graphene is controled by the 
backgate with a similar principle. The most interesting fea- 
ture here is that the locally gated region can be controled 
independently due to the screening of the embedded local 
gate, as is evident in both Figs. 5(b) and 5(c). This indepen- 
dent control leads to the four quadrants of the conductance 
map with the two boundaries perpendicular to each other 
[13], as contrary to those observed in top-gated devices [4, 
5,6,7,8,11,12]. 



4.3 Graphene superlattices 

Next we turn to the possibility of generating a graphene su- 
perlattice by fabricating a series of patterned topgates. As 
sketched in Fig. 6(a), the PDE problem is defined within 
the region above the graphene sheet in order to solve the 
electric potential due to the topgates with various voltages 
y tg = —5, —3, • • • ,5 V. The contribution from the backgate 
is assumed to be uniform and can be treated independently. 
The strategy here is to compute first the carrier densities due 
to topgates, and then include the backgate contribution to 
yield the total carrier density that finally gives the energy 
band offset profile from Eq. (28). In this case the approxima- 
tion (27) is rather acceptable, taking the voltage difference 
simply as V tg fixed as position-independent. Hence we will 
compare the results from all of the three approaches. 

The computed carrier densities, «pd by PDM, tiqq by 
QCM, and «c the classical contribution, are shown in Fig. 
6(b). The curves of «pd and «qc almost coincide with each 
other. The relatively thick 40 nm of AI2O3 suppresses the 
quantum correction to a reasonably small amount, such that 
here the CCM is not a bad approximation, either. The dis- 
crepancy between «c and «pd (or «qc) is less pronounced at 
the regions between each adjacent pair of topgates since the 
quantum correction \An\ roughly decreases with the 3 /2-th 
power of the distance to the gate, as mentioned in Eq. (26). 

The carrier density modulation follows the patterned 
topgates with a periodicity of lOOnm, giving rise to a pe- 
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Fig. 6 Formation of a graphene superlattice using a series of patterned 
topgates and a global backgate, which are treated separately, (a) It- 
erative solution u(x,z) with V tg = 5 V, considering only the topgates. 
(b) The corresponding carrier densities based on all of the three ap- 
proaches, (c) The profiles of the energy band offset from the total car- 
rier density composed of the patterned topgates and the uniform back- 
gate contributions. 



riodic potential profile V(x) as shown in Fig. 6(c), where a 
backgate contribution with V^g = 30 V is taken into account. 
Since V(x) is related to n(x) through a square-root relation, 
Eq. (28), the shape of V(x) can be a bit different from that 
of n(x), which is similar to a sine-like wave, especially for 
those V (x) that alternate between positive and negative val- 
ues. Note that the backgate voltage chosen in Fig. 6(c) re- 
sults in a rather symmetric V tg = — 3 V curve since the cor- 
responding carrier density n(x) alternates symmetrically be- 
tween positive and negative. In general, the alternation of 
n(x) is not necessarily symmetric (about the charge neutral- 
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ity point n = 0), and the resulting V(x) profile can be of 
peculiar shapes. 



4.4 Linear potential 

In the previous example, we have demonstrated that fabri- 
cating a series of patterned topgates may generate a periodic 
potential, which, combined with a potential linear in posi- 
tion, forms the prerequisite of observing the Bloch-Zener 
oscillation in graphene [22]. In the last demonstrating ex- 
ample here, we point out a simple way to generate the linear 
potential: using a tilted backgate. As sketched in Fig. 7(a), 
where we consider a position-varying thickness of Si02 with 
a slope of s = 0.05 (an increase of 50 nm per micron). In this 
case the quantum correction does not play a role, and we 
show in Fig. 7(b) only the carrier densities from the CCM 
and the PDM, which coincide to each other. 

Since the classical capacitance model works well here, 
with the oxide thickness d(x) = do + sx, where do = 300nm 
is the thickness at the center, we can describe the carrier 
density as n(x) = eV g /e(do + sx) and hence the potential as 

V(x) = — sgn(V g )h\'F ne \V g \ /e(do + sx). The slope of the 
potential at x = 0, dV (x) / dx\ x=0 , together with the intercept 
V(x = 0), allows us to approximate the potential with a lin- 
ear model, 



V(x)xsV + Sx 
V = -sgn(V g ) 

S = sgn(V g ) 



graphene 




x 0.274 25 eV 



(29) 



.3/2 



x 0.137 13eVnm 



-l 



In Fig. 7(c), we plot the potential profiles obtained from «pd 
and from the linear model given by Eq. (29); the consistency 
is almost perfect. 



5 Summary 
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Fig. 7 (a) Electrostatic potential u(x,y) inside a trapezoidal oxide 
layer. The tilted backgate with a slope of 0.05 generates a carrier den- 
sity (b) that varies almost linearly with position. The corresponding 
potential profiles (c) also exhibit a linear behavior. 



the usage of the Matlab pdetool, the major part of this work 
provides a self-contained instruction to calculating the car- 
rier density of bulk graphene subject to any kind of gating 
geometry. 

To demonstrate the application of the gate-induced car- 
rier density modulation in graphene, practical examples 
have been illustrated, including the graphene pnp junction 
by using an embedded local gate in additional a global back- 
gate [13], graphene superlattice potential by a series of pat- 
terned topgates, and the quasi-linear potential by using a 
tilted backgate. These examples correspond to the experi- 
mental conditions that provide a flexible platform to test the 
physics of Klein backscattering [20, 12,21] and the Bloch- 
Zener oscillation [22] in graphene. 
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In conclusion, theories of the gate-induced carrier density 
modulation in bulk graphene have been reviewed, including 
the classical capacitance model that most of the time plays 
the dominant role, the Poisson-Dirac method that takes into 
account the quantum correction by paying the price of iter- 
ation, and the quantum capacitance model that can analyti- 
cally deal with the quantum correction. The analytical QCM 
and the numerically iterative PDM are shown to be exactly 
equivalent to each other for the uniform capacitor case. For 
the case of nonuniform capacitors, the QCM is also shown 
to work equally well with PDM, as long as the approxima- 
tion (27) is not ill-defined. Along with the brief summary for 



References 

1 . A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and 
A. K. Geim, "The electronic properties of graphene," Rev. Mod. 
Phys., vol. 81, p. 109, 2009. 

2. S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, "Electronic 
transport in two-dimensional graphene," Rev. Mod. Phys., vol. 83, 
pp. 407^170, May 2011. 

3. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, 
S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, "Electric field 
effect in atomically thin carbon films," Science, vol. 306, no. 5696, 
pp. 666-669, 2004. 

4. B. Huard, J. A. Sulpizio, N. Stander, K. Todd, B. Yang, and 
D. Goldhaber-Gordon, "Transport measurements across a tunable 



12 



potential barrier in graphene," Phys. Rev. Lett., vol. 98, p. 236803, 
Jun 2007. 

5. J. R. Williams, L. DiCarlo, and C. M. Marcus, "Quantum hall 
effect in a gate-controlled p-n junction of graphene," Science, vol. 
317, no. 5838, pp. 638-641, 2007. 

6. B. Ozyilmaz, P. Jarillo-Herrero, D. Efetov, D. A. Abanin, L. S. 
Levitov, and P. Kim, "Electronic transport and quantum hall effect 
in bipolar graphene p-n-p junctions," Phys. Rev. Lett., vol. 99, p. 
166804, Oct 2007. 

7. G. Liu, J. J. Velasco, W. Bao, and C. N. Lau, "Fabrication of 
graphene p-n-p junctions with contactless top gates," Appl. Phys. 
Lett, vol. 92, no. 20, p. 203103, 2008. 

8. R. V. Gorbachev, A. S. Mayorov, A. K. Savchenko, D. W. Horsell, 
and F. Guinea, "Conductance of p-n-p graphene structures with 
air-bridge top gates," Nana Letters, vol. 8, no. 7, pp. 1995-1999, 
2008, pMID: 18543979. 

9. V. V. Cheianov and V. I. Fal'ko, "Selective transmission of 
Dirac electrons and ballistic magnetoresistance of n-p junctions 
in graphene," Phys. Rev. B, vol. 74, no. 4, p. 041403, JUL 2006. 

10. M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, "Chiral tun- 
nelling and the klein paradox in graphene," Nat. Phys., vol. 2, 
no. 9, p. 620, SEP 2006. 

11. N. Stander, B. Huard, and D. Goldhaber-Gordon, "Evidence for 
klein tunneling in graphene p-n junctions," Phys. Rev. Lett., vol. 
102, p. 026807, Jan 2009. 

12. A. F. Young and P. Kim, "Quantum interference and klein tun- 
nelling in graphene heterojunctions," Nat. Phys., vol. 5, no. 3, pp. 
222-226, MAR 2009. 

13. S.-G. Nam, D.-K. Ki, J. W. Park, Y. Kim, J. S. Kim, and H.-J. 
Lee, "Ballistic transport of graphene pnp junctions with embedded 
local gates," Nanotechnology, vol. 22, no. 41, p. 415203, 2011. 

14. J. Guo, Y. Yoon, and Y. Ouyang, "Gate electrostatics and quantum 
capacitance of graphene nanoribbons," Nana Letters, vol. 7, no. 7, 
pp. 1935-1940, 2007. 

15. J. Fernandez-Rossier, J. J. Palacios, and L. Brey, "Electronic 
structure of gated graphene and graphene ribbons," Phys. Rev. B, 
vol. 75, p. 205441, May 2007. 

16. T. Fang, A. Konar, H. Xing, and D. Jena, "Carrier statistics and 
quantum capacitance of graphene sheets and ribbons," Applied 
Physics Letters, vol. 91, no. 9, p. 092109, 2007. 

17. A. A. Shylau, J. W. Klos, and I. V. Zozoulenko, "Capacitance 
of graphene nanoribbons," Phys. Rev. B, vol. 80, p. 205402, Nov 
2009. 

18. T. Andrijauskas, A. A. Shylau, and I. V. Zozoulenko, "Thomas- 
fermi and poisson modeling of gate electrostatics in graphene 
nanoribbon," Lithuanian Journal of Physics, vol. 52, no. 1, pp. 
63-69, 2012. 

19. S. Luryi, "Quantum capacitance devices," Applied Physics 
Letters, vol. 52, no. 6, pp. 501-503, 1988. 

20. A. V. Shytov, M. S. Rudner, and L. S. Levitov, "Klein backscat- 
tering and fabry-perot interference in graphene heterojunctions," 
Phys. Rev. Lett., vol. 101, p. 156804, Oct 2008. 

21. M.-H. Liu and K. Richter, "Efficient quantum transport simulation 
for bulk graphene heterojunctions," Phys. Rev. B, vol. 86, p. 
115455, Sep 2012. 

22. V. Krueckl and K. Richter, "Bloch-zener oscillations in graphene 
and topological insulators," Phys. Rev. B, vol. 85, p. 115433, Mar 
2012. 

23. Partial Differential Equation Toolbox™ User's Guide, Matlab 
2012a ed., The MathWorks, Inc., 2012. 

24. E. I. Rashba, "Properties of semiconductors with an extremum 
loop i. cyclotron and combinational resonance in a magnetic field 
perpendicular to the plane of the loop," Sow Phys. Solid State, 
vol. 2, p. 1109, 1960. 

25. Y. A. Bychkov and E. I. Rashba, "Properties of a 2d electron-gas 
with lifted spectral degeneracy," JETP Lett., vol. 39, p. 78, 1984. 



Ming-HaoLiu (IPJ*) 



26. M. Gmitra, S. Konschuh, C. Ertler, C. Ambrosch-Draxl, and 
J. Fabian, "Band-structure topologies of graphene: Spin-orbit 
coupling effects from first principles," Phys. Rev. B, vol. 80, p. 
235431, Dec 2009. 

27. S. Abdelouahed, A. Ernst, J. Henk, I. V. Maznichenko, and 
I. Mertig, "Spin-split electronic states in graphene: Effects due to 
lattice deformation, rashba effect, and adatoms by first principles," 
Phys. Rev. B, vol. 82, p. 125424, Sep 2010. 

28. A. Yamakage, K.-I. Imura, J. Cayssol, and Y. Kuramoto, 
"Interfacial charge and spin transport in Fi topological 
insulators," Phys. Rev. B, vol. 83, p. 125401, Mar 2011. 

29. H. Y. Tian, Y. H. Yang, and J. Wang, "Interfacial charge current 
in a magnetised/normal graphene junction," European Physical 
Journal B, vol. 85, no. 8, AUG 2012. 

30. A. Yamakage, K. I. Imura, J. Cayssol, and Y. Kuramoto, "Spin- 
orbit effects in a graphene bipolar pn junction," EPL, vol. 87, no. 4, 
AUG 2009. 

31. M.-H. Liu, J. Bundesmann, and K. Richter, "Spin-dependent klein 
tunneling in graphene: Role of rashba spin-orbit coupling," Phys. 
Rev. B, vol. 85, p. 085406, Feb 2012. 

32. M. Rataj and J. Barnas, "Graphene p-n junctions with nonuniform 
rashba spin-orbit coupling," Appl. Phys. Lett., vol. 99, no. 16, p. 
162107, 2011. 

33. J. Sanchez-Barriga, A. Varykhalov, M. R. Scholz, O. Rader, 
D. Marchenko, A. Rybkin, A. M. Shikin, and E. Vescovo, "Chem- 
ical vapour deposition of graphene on Ni(l 1 1) and Co(0001) and 
intercalation with Au to study Dirac-cone formation and Rashba 
splitting," Diamond and Related Materials, vol. 19, no. 7-9, SI, 
pp. 734-741, JUL-SEP 2010, 20th European Conference on Dia- 
mond, Diamond-Like Materials, Carbon Nanotubes and Nitrides, 
Athens, GREECE, SEP 06-10, 2009. 

34. J. Xia, F. Chen, J. Li, and N. Tao, "Measurement of the quantum 
capacitance of graphene," Nature Nanotechnology, vol. 4, no. 8, 
pp. 505-509, AUG 2009. 

35. S. Droscher, P. Roulleau, F. Molitor, P. Students, C. Stampfer, 
K. Ensslin, and T. Ihn, "Quantum capacitance and density of 
states of graphene," Applied Physics Letters, vol. 96, no. 15, p. 
152104, 2010. 

36. L. A. Ponomarenko, R. Yang, R. V. Gorbachev, P. Blake, 
A. S. Mayorov, K. S. Novoselov, M. I. Katsnelson, and A. K. 
Geim, "Density of states and zero landau level probed through 
capacitance of graphene," Phys. Rev. Lett., vol. 105, p. 136801, 
Sep 2010. 

37. J. Martin, N. Akerman, G. Ulbricht, T. Lohmann, J. H. Smet, 
K. Von Klitzing, and A. Yacoby, "Observation of electron-hole 
puddles in graphene using a scanning single-electron transistor," 
Nature Physics, vol. 4, no. 2, pp. 144-148, FEB 2008. 

38. H. Xu, Z. Zhang, and L.-M. Peng, "Measurements and 
microscopic model of quantum capacitance in graphene," Applied 
Physics Letters, vol. 98, no. 13, p. 133122, 2011. 



